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ABSTRACT 

The merger and accretion probabilities of dark matter halos have so far only been calculated for an 
infinitesimal time interval. This means that a Monte-Carlo simulation with very small time steps is 
necessary to find the merger history of a parent halo. In this paper we use the random walk formalism 
to find the merger and accretion probabilities of halos for a finite time interval. Specifically, we find 
the number density of halos at an early redshift that will become part of a halo with a specified final 
mass at a later redshift, given that they underwent n major mergers, n = 0,1,2,... . We reduce 
the problem into an integral equation which we then solve numerically. To ensure the consistency of 
our formalism we compare the results with Monte-Carlo simulations and find very good agreement. 
Though we have done our calculation assuming a flat barrier, the more general case can easily be 
handled using our method. This derivation of finite time merger and accretion probabilities can be 
used to make more efficient merger trees or implemented directly into analytical models of structure 
formation and evolution. 

Subject headings: cosmology:theory-dark matter 



1. INTRODUCTION 

Since the introduction of the spherical collapse model 
of dark matter halos by Press and Schcchtcr (Press & 
Schechter 1974) and its generalizations such as extended 
Press-Schechter (EPS) (Bond et al. 1991) and ellipsoidal 
collapse (Sheth & Tormen 2002) it has been used ex- 
tensively in the cosmological literature as a fast and ac- 
curate method to quantify the distribution of collapsed 
dark matter objects in the universe. This in turn is the 
backbone of semi-analytic theories of galaxy formation 
and evolution (Cole et al. 2000). Modified versions of 
the excursion set formalism underlying EPS have also 
found application in different contexts such as ionized 
bubble growth in the early universe (Furlanetto et al. 
2004) . The advantages of this method compared to di- 
rect N-body simulation include its superior speed, which 
allows the exploration of large ranges of parameter space, 
and redshift range than is currently accessible to N-body 
simulations. 

In the EPS formalism one finds the probability 
P{M\, zi\M2, Z2)dMi which gives the probability of a 
point mass being part of a halo with mass in between 
Mi and M% + dM\ at redshift Z\ given that it was (or 
will be) part of a halo with mass at redshift z%. Dur- 
ing this redshift interval it could merge with any num- 
ber of halos whose masses add up to \M\ — Ma|. How- 
ever, from the viewpoint of galaxy formation, accretion 
of small halos into the larger one do not have the ability 
to change the evolution of the galaxy or galaxies inside 
it. Indeed, it is assumed that only at major mergers in 
which the mass of the merged halo satisfies a condition 
to be large enough, have this ability. It is therefore in- 
teresting to ask the question: given the criterion for a 
major merger, is it possible to find an analytical result 



to describe the progenitor distribution of a parent halo, 
based on how many times they have undergone a ma- 
jor merger? For example, what is the number density of 
halos at redshift, say, z — 1 which underwent n major 
mergers, n = 0, 1, 2, and eventually ended up being a 
galaxy halo at the present time. 

Here we present an analytical method based exclusively 
on excursion set assumptions to find these number den- 
sities. We reduce the problem to an integral equation 
which can then be solved numerically. Our derivation is 
given in section 2. We show in section 3 that Monte- 
Carlo (MC) simulations agree with the results of our 
semi-analytic model. 

It is also possible to cast the integral equation into 
a scale invariant form. This can be a great advantage 
since we just need to solve the integral equation once 
and scale the result to find the general formula. This is 
done in section 4. We discuss how the method presented 
in this paper can be generalized and conclude in section 
5. 

A few appendixes are added for further clarification. 
Appendix A gives a brief introduction to the excursion 
set formalism and defines our notation. In Appendix 
B we give a straightforward method to solve the integral 
equation. Finally in Appendix C we describe very briefly 
how we implemented our Monte-Carlo simulation. 

2. THE ACCRETION PROBABILITY 

The main task of this paper is to find an analytical 
result for f a cc(S2\^2, Si, u>±), the probability that a ran- 
dom walk starting from (Si,wi) has its first upcrossing 
between 5*2 and 5*2 + dS*2 at the barrier height 0J2 given 
that it never had a jump larger than M res between wi 
and u>2- Here S — a 2 (M) is the rms mass fluctuation 
inside spheres of mass M , and w is a monotonically in- 



ai+ dw 




Fig. 1. — This figure shows a sample trajectory for a point that 



belongs to a halo of mass Sf at time u>f in a S vs 



diagram. 



From this figure we can find places where the halo to which this 
point belongs splits into two smaller halos of masses larger than 
M reB by finding the places where the first upcrossing of the random 
walk has a jump larger than AS = S(M) — S(M — M res ) when 
the barrier height is increased from utow + dw. In this figure the 
mass resolution is such that there is just one jump large enough 
to be considered a merger, and occurs where the first upcrossing 
jumps from Si to S2. The rest of the time the halo will just 
accrete masses smaller than M rea . We want to find the fraction 
of random walks that start from (Sf,Wf) and will have their first 
upcrossing between Si and Si + dSi at the barrier height u>i and 
will have one and only one large enough jump to be considered 
as merger. We first find the fraction of the random walks that, 
starting from (Sf,u>f), just accrete mass and end up somewhere 
in the mass range (Si, Si + dSi) at the barrier height u>. Then 
we take the random walks that passed the previous test and find 
the fraction of them that have a merger from Si to some S2 > 
S(M(Si) — M rea ) during an infinitesimal time interval duj. Finally, 
from these random walks we take the fraction that, starting from 
(S2,w), only accrete and end up in the mass range (Si, Si + dSi) 
at the barrier height Wi. The outcome of this is the equation 2. 
To find the total fraction of the random walks that have one and 
only one merger we need to integrate over intermediate values of 
Si, S2 and w, keeping in mind that S2 cannot come closer to Si 
than what the mass resolution allows us. This give us the equation 
3. 



creasing function of redshift which for the case of an EdS 
universe takes the familiar form of 1.68 x (1 + z) ' . 
In order to accomplish this we let a number of random 
walks start from (S/,w/) (see fig. 1), where in this pa- 
per Sf and ujf denote a 2 (Mf) and co(zf) respectively, 
and so forth for other subscripts. Setting a barrier at 
Ui > ujf, that is at earlier times, we know the fraction 
of random walks that have their first upcrossing in the 
interval (Si, Si + dSi), regardless of whether they accrete 
or merge, is (see equation A4): 



ftot{Si\u>i, Sf,uj f ) dS. : 

(uii - Ulf) 
: (2n)V*(Si-S f )V* 



exp 



(u>j - LOf) 2 

'2(Si-S f ) 



dSi (1) 



The next step is to notice that the fraction of random 
walks that start from (S/,w/) and have their first up- 
crossing between Si and Si + dSi at the height u>i and 
have one and only one jump from Si to S2 during the 

' For further explanation of excursion set formalism, see Ap- 
pendix A 



interval u> and u> + dw is: 

focc(Si\w,Sf,Wf)dSi x f(Si -> S 2 ;uj)dS 2 du) 

X facc(Si\Ui,S 2 ,u)dSi, (2) 

where /(Si — > S 2 ;u>) dS 2 dw is the probability that a ran- 
dom walk will have a sudden jump from S\ to somewhere 
between S 2 and S 2 + dS 2 in an infinitesimal interval duj. 
Its form for the Spherical collapse model is given by equa- 
tion A5. 

Now, if we integrate over all Si, S 2 and u) keeping in 
mind that M(Si) — M(S 2 ) must be larger than M res to 
be classified as a merger, we find the fraction of walks 
that have undergone one and only one merger to be: 

rf>i rSi-A'(Si) rSi 

flmerger(Si\&i,Sf,LOf)= / dhJ dS\ I 

Juj f JSf JSi+A(Si) 

f acc (Si\uj,Sf,ujf)f(Si -> S 2 ;uj)facc(S i \uj ll S 2 ,u) 

Here A (Si) = cr 2 (M(Si)-M res )-Si is the closest dis- 
tance that S 2 can be brought to Si and still have a jump 
large enough to be considered as a merger. Analogously 
define A' (Si) = Si - a 2 (M (Si) + M res ). 

We can similarly find the fraction of random walks 
that have undergone two and only two mergers, 
f 2 mergers(Si\u>i, Sf,C0f). We get an expression like equa- 
tion 3 but with six integrations, and higher terms involve 
more integrations yet. Adding all of these terms must 
give: 

f tot (Si\wi,Sf,Uf) = f acc (Si\uji,Sf,ujf) + 

flmerger(Si\0Ji, Sf,U)f) + f 2m ergers(Si\oJi, Sf,U)f) H 

(4) 

which says that the fraction of random walks starting 
from (Sf,ujf) and passing through a barrier at cj, for 
the first time between (Si, Si + dSi), no matter what 
happened during their journey, is equal to the sum of 
the fractions of walks which underwent no mergers, one 
merger, two mergers, etc. 

Putting the integral formulae for f\merge(Si\uJi, Sf,u>f) 
and higher order terms into equation 4, we will find 
an integral equation for the unknown f a cc(Si\uJi, Sf, cjf). 
However, this equation written in this form is not compu- 
tationally tractable since higher order terms will require 
a prohibitive number of integrations. 

Note that the conditional probability densities in 
Eq. (4) are analogous to propagators. It is instructive 
to visualize each term in equation 4 using a diagram- 
matic notation, as in Fig. 2. The figure makes it clear 
that this equation expands the full propagator in terms of 
bare propagators (the accretion probability) with inter- 
actions (mergers). We can use this insight to rearrange 
the above equation as in the figure, which shows that one 
can resum the terms in Eq. (4) to write down a tractable 
integral equation. 

Writing the second equality of figure 2 in the language 
of equation 4, we find the important result: 

ftot(Si\iOi, Sf, LOf) = facc(Si\uJi, Sf,UJf) 
■Mi rS t -A'(Si) j-Si 

dui / g?Si / dS 2 

f JS f /Si+A(Si) 

f acc (S 1 \uj,Sf,ujf)f(S 1 -> S 2 ;uj)f tot (S t \uji,S 2 ,uj)(5) 



X- 



-x---x--- + 



Fig. 2. — Feynman diagrams illustrating equation 4 for the total 
probability for a random walk starting from (Sf,Wf) to have its 
first upcrossing between Si and Si + dSi at u>i . A solid line denotes 
a period of time in which any sequence of accretion and merger 
events can take place. A dashed line denotes a period of time in 
which only accretion takes place and a cross indicates a merger. 
The first equality says that the total probability is equal to the 
sum of the fractions of walks which underwent no mergers, one 
merger, two mergers, etc. We can rearrange the sum to give the 
second equality, which states that the total probability is equal to 
the sum of the probability to have no mergers and the probability 
to have at least one merger. 
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Fig. 3. — This figure shows p a cc versus Si, with parent mass 



M 



f 



M 



f 



and redshift 



2 X fO 12 Mo, mass resolution M rea — -jg- 
Zf = 0. The panels show ptot (solid line), our analytic result 
for pace (dashed line), and the MC result for p aC c (histograms) 
for increasing redshift from top-left to bottom-right: z = 0.184, 
z = 0.267, z = 0.348 and z = 0.5. 



This is a Voltera integral equation which we will solve nu- 
merically. We refer the interested reader for a discussion 
of solving this equation numerically to appendix B. In 
the next section we compare these semi-analytic results 
to Monte-Carlo simulations. 

3. COMPARISON WITH MONTE-CARLO SIMULATION 

There are several methods for generating a merger tree 
using Monte-Carlo techniques (see Somerville & Kolatt 
(1997) and Cole et al. (2000)). Here we choose to use 
a binary merger tree with accretion method mainly due 
to its simplicity of implementation. In this scheme the 
time interval is chosen to be so small that the probability 
of a merger is very low. This in turn ensures that the 
probability of more than one merger is negligible in a 
given time step. Then a halo one time step back will 
have a smaller mass due to accretion or division into two 
progenitors. The details of the method can be found in 
Appendix C. 

It should be noted that our Monte-Carlo (MC) sim- 
ulation gives number weighted probabilities. However, 
face in formula 5 is a mass weighted probability. We can 
easily change f aC c to a number weighted probability p acc 

Mi r 

USmg p acc = -jfcfacc. 

Each panel of figure 3 shows p acc versus Si for a dif- 
ferent lookback redshift zi, with parent mass Mt = 



Fig. 4. — The dashed lines in these figures show our analytic 
result for p acc vs. Si with the same Mf and Zf as figure 3, but 
here we look at the distribution at a fixed lookback redshift Zi = 0.5 



for decreasing Mres' from top-left to bottom-right Mre 



M 



M r 



M re 



and Mr- 



Solid lines show ptot 



and histograms the results of MC simulations for p a c 



2 x 10 12 Af Q and redshift Zt = 0. The mass resolution 

is fixed at M res = -j^-. The dashed line is our analyt- 
ical solution, the histogram shows the result of the MC 
simulation and the solid line is ftot given by equation 1. 
We see an excellent agreement between the MC simula- 
tion and our analytic result, and notice some intuitively 
sensible trends in the figures that are worth mentioning. 
First, p acc — ptot for Si < 4.1. That is because to have 
a merger we need Mi < Mf — M res = 1.8 x 10 12 , which 
corresponds to Si — 4.1. For Sj smaller than this, the 
mass jump is always less than mass resolution and only 
accretion can happen; hence, p aC c = Ptot in this region. 
For Si larger than 4.1 mergers are allowed, so the proba- 
bility of having one or more merger is non-zero and p aC c 
will be less than p tot . Also, the probability of having 
at least one merger, p to t — Pace, increases monotonically 
with increasing redshift as more and more halos have a 
chance to undergo a merger. Finally, as we look further 
back in time, halos with smaller masses have a chance 
to reach Mf by just accreting so p aC c spreads to smaller 
masses with increasing redshift. 

In figure 4 we show p acc vs. Si with the same Mf 
and zj as above, but here we look at the distribution 
with different choices of M res at a fixed lookback redshift 
Zi = 0.5. Again the agreement between our analytical re- 
sult and the MC simulations is very good. Again, there 
are some intuitively reasonable trends in the figure that 
should be noted. As we discussed above, p aC c must be 



equal to p to t for M, > M\ 



limit 



Mf - M res . As M r 



gets smaller, this limiting mass becomes larger. There- 
fore Surmt = S(M Hmi t), the Si below which p acc = p tot , 
becomes smaller, which can easily be seen in the figure. 
Notice that different panels have different scales. Also 
given a fixed lookback redshift, decreasing the mass reso- 
lution increases the number of events we classify as merg- 
ers thus raising the probability for a halo to have at least 
one merger. Since ptot is not affected by the choice of 
M res , Pacc accordingly decreases with decreasing M res - 



4. SCALED SOLUTION 

In general, for any given prescription for M res one 
can find the solution of integral equation 5 to obtain 
facc{Si\^i, Sf,ujf). Generally, this needs to be solved for 
each given final halo mass Mf. However, if we impose a 
special mass resolution for a chosen cosmology it is pos- 
sible to cast the integral equation into a scale invariant 
form. Then we need only solve this equation once. To 
achieve this, we need to define M res so as to satisfy the 
following two equations simultaneously: 



A(S) = a 2 (M(S) - M res ) -S = CxS 



(6) 



and 



A'(S) = S-a 2 {M{S) + M res ) = C x S (7) 



where C and C are constants independent of S. Recall 
that A(S) and A'(S') appear in the limits of integration 
of equation 5. For M res small compared to M(S) we 
can easily see, by Taylor expansion, that to second order 
in M res /M these equations can be satisfied if we take 
C = C and the mass resolution as 



M res {S) = - 



dS/dM 



C. 



(8) 



For example for a scale invariant matter power spectrum 
with power index n, P(k) oc k n , the above equation gives 



M res = 



-CM, 



parent 



(9) 



i.e. a merger is defined when the mass of any progenitor 
of the halo is larger than a constant fraction of the parent 
mass 1 . 

With this criterion for M res we are ready to rewrite 
the integral equation 5 in a scale invariant form. To do 
so we define the new variables 



u ee S/S f 

e = (u-uj f )/sy 2 

and the functions f acc and / acc : 

facc(Sl\u,Sf,U}f) = S J 1 , f acc ^-\ 



(10) 



facc{ui\6) 



{2n) 1 / 2 (u 1 - l) 3 / 2 



exp 



2(ui - 1) 



_ ( 27 r)V2 K - 1)3/2 

a q2 



(2 7 r) 1 /2( Ul _ 1)3/2 eXP 2 ( Ul -l) 

1 



(2tt) 1 /2( U2 _ Ul )3/2 (27r)V2( Ui _ U2 )3/2 * 2 { Ul - Ui) 

(13) 

For a given C this equation can be solved for f a cc(u\9). 
This calculation can be facilitated by noticing that the 
integral over S 2 can be done analytically. Having found 
facc{u\0) one can find f acc (Si\wi, Sf,u> f ) for arbitrary 
Si, uji, Sf and u)f using equations 11. 

The result of this calculation for a power law matter 
power spectrum with n = —1 and C = 0.01 is shown 
in figure 5. The solid line, as usual, denotes ptot, the 
lighter histogram in each panel indicates the result of 
MC simulation for p acc and the dashed line on top is 
our numerical solution of equation 12 scaled according 
to equation 11. The panels from top- left to bottom right 
are for ui — u)f — 0.2, 0.3, 0.4 and 0.5. Assuming u> — 
1.69(1 + z) and z; = these correspond to lookback 
rcdshifts 0.118, 0.178, 0.237 and 0.296 respectively. One 
can see that the agreement is very good for all rcdshifts 
considered. 

The darker histogram in figure 5 shows the result of 
our MC simulation for pimerger(Si, u>i, Sf, u>f). This is 
the number density of halos in a given range (Si, Si+dSi) 
that have a parent halo of mass corresponding to Sf at 
time u)f and have merged once in their journey from uoi 
to uif. Now that we have found f acc it is possible to cal- 
culate Plmerger{Si, U>i,Sf,U)f) USmg equation 3. Plmerger 



is nothing but /; 



Imerger 



multiplied by ^f- to convert from 



With these definitions equation 5 can be written in the 
manifestly scale invariant form 



r@i rUi(l — C) rut 
1 = facc(Ui\0i) + d6 dui 

JO Jl Ju!(l 



r u t (l-C) 

dui I du 2 

i(l+C) 

\0)K{6i,Ui,ui,U2,6) (12 
where the kernel for the spherical collapse model is 



mass density to number density. This result is shown by 
the dot-dashed line on top of the histogram. The match 
is very good. As expected, for small Az the probability 
of one merger is much smaller than probability of accre- 
tion, which can be seen in the top-left plot. Also in this 
plot, we can see that the tail of the distribution pimerger 
approaches p tot . This says that the probability of hav- 
ing more than one merger is negligible for small Az, as 
expected. 

On the other hand, when Az gets larger more and 
more halos have a chance to merge, so p acc flattens and 
Pimerger rises. Also, with a large Az there is a finite 
facc{ui\0) probability of having more than one merger since the tail 
of Pimerger is considerably below p to t ■ One can continue 
(ll)this calculation and find P2mergers, Pzmergers and so on, 
which we have not shown here. Where this hierarchy 
should be terminated clearly depends on how far we look 
back in time: a larger Az means more chance of a merger 
and therefore requires higher merger terms. 



1 Note that more general solutions exist if C 7^ C which can pos- 
sibly be used to loosen the relationship between M res and M par ent- 
It seems physically reasonable to choose M re s much smaller than 
and proportional to M paren t which is the case when C = C . 



5. DISCUSSION AND CONCLUSION 

We have used the random walk formalism to find the 
accretion probability, i.e. the probability for a parent 
halo to have a progenitor in a given mass interval at 
a given earlier time given that it has not merged with 
a halo of a mass larger than the mass resolution. As 
a concrete example we have worked out the accretion 
probability in the special case where the barrier is flat, 
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Fig. 5. — Halo accretion and merger probabilities for a power 
law matter power spectrum with n = — 1 are shown here. We take 
C = 0.01 where C is denned in equation 9. The solid line denotes 
Ptot, the lighter histogram in each panel indicates the result of MC 
simulation for p a cc and the dashed line is our numerical solution of 
equation 12 scaled according to equation 11. The darker histogram 
shows the result of our MC simulation for pi m erger(Si,LOi, Sf,ujf), 
the number density of halos in a given range (Si, Si + dSi) that 
have a parent halo of mass corresponding to Sf at time Uf and 
have merged once in their journey from u>i to uif. The panels 
from top-left to bottom right are for u> — u)j = 0.2, 0.3, 0.4 and 
0.5. Assuming u> = 1.69(1 + z) and zj = these correspond to 
lookback redshifts 0.118, 0.178, 0.237 and 0.296 respectively. 



the mass resolution is constant and we look backward 
in time. However this method can be extended to solve 
more general problems. 

For example, while we have used a constant barrier, 
it is well known that this barrier shape does not match 
the results of N-body simulations. However, our formal- 
ism can be generalized to the case of a moving barrier, 
which has proven to give a very good match to N-body 
simulations. One only needs to find the appropriate for- 
mulae for ftot and f(S\ — > Sai to) for the moving barrier 
and solve the integral equation 5. These functions can 
in general be found numerically, using for example the 
method of Zhang & Hui (2006). However there are ana- 
lytical results for simple barriers that reproduce the re- 
sults of N-body simulations, e.g. the square root barrier 



(Mahmood & Rajesh 2005; Giocoli et al. 2007). Since 
the aim of this paper is not to compare with numerical 
simulations, we will leave this calculation for future work. 

Here we have always looked backward in time. How- 
ever, in certain cases it might be more convenient to find 
f acc in the forward sense. In that case it gives the prob- 
ability of a halo at an early time being accreted by a 
larger halo in a given mass interval. This problem can 
be solved with our formalism by using the forward form 
of equation 5, with the forward form of / tot and f merge 
(Lacey & Cole 1993). Given a population of objects of 
mass M at time t, some of these objects will be destroyed 
in the course of their evolution by merging with other ob- 
jects. To find the fraction of objects that have survived 
from the initial time to the observation time (see Verde 
et al. 2001, for the case of clusters of galaxies) we need to 
find the fraction of objects whose halos have not merged 
with halos more massive than a given threshold; in other 
words they have only accreted from the initial redshift 
to the redshift of observation. This is precisely what our 
formalism calculates. 

Finally, this method can lead to a major improvement 
in the speed of merger tree generation in Monte-Carlo 
simulations. Since there was no formula for the accre- 
tion or merger probabilities in a finite time interval, past 
MC codes had to use infinitesimal time-steps to be able 
to use the known formula for merger probabilities (eqn. 
A5), making the computation very time consuming. In 
this paper we have presented methods to calculate the 
accretion and merger probabilities for any given time in- 
terval which can be used to generate trees with larger 
time-steps and hence in less computational time. 
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APPENDIX 

BRIEF OVERVIEW OF HALO MODEL 

In the excursion set formalism one assumes all the matter in the universe is inside collapsed objects, halos, and the 
aim is to find the number density of these objects for different halo masses. In order to do that one starts from the 
initial density field of dark matter, which is assumed to be Gaussian, and finds its present time distribution using linear 
theory. The next step is to take non- linear effects into account, but non-linear theories are generally very complicated 
and hard to calculate analytically. The halo model circumvents this difficulty by usage of the simplest non-linear 
model, i.e. the spherical collapse (SC) model, for which a simple analytical solution exists. In the SC model one finds 
that if an overdense sphere collapses at redshift z then its initial overdensity extrapolated linearly to the present time 




Fig. A6. — This graph shows the rms mass fluctuation inside spheres of mass M, i.e S(M), vs mass for our fiducial ACDM cosmology. 



will have the value S c D(0)/ D(z) where 5 C depends weakly on cosmology and D is the linear growth factor. This fact 
is used in the model by taking a sphere centered on any point in space and finding these two quantities inside this 
sphere: 

S = tzl (Al) 



and 



S 



a 2 oc 



k 2 P{k)W 2 (kR)dk 



(A2) 



where p is the mean density of the universe and p the mean density inside the sphere. Notice that all of the quantities 
are calculated using the initial density linearly extrapolated to the present time. P is the matter power spectrum, 
W is the top-hat window function in real space and R = (j^) 1 . The constant of proportionality is found from er g . 
For hierarchical cosmologies S is a decreasing function of mass inside the sphere and goes to zero for large radii or 
equivalently for large masses. Therefore, there is a one to one map between S and M. In figure A6, we give a plot of 
S as a function of M for the ACDM cosmology with fi m = 0.25, h = 0.73 and as — 0.9. 

The property that makes this formalism so appealing is that in the case of a Gaussian initial distribution by reducing 
the radius of the sphere, S will execute an uncorrelated random walk for which S is the time like quantity. The plot 
for S versus S for a random point in space will be a realization of a 1-D random walk (see figure 1). To find the mass 
of the halo to which this point belongs at redshift z, one draws a horizontal line at S c (z) and finds the point which 
cuts the random walk for the first time. S at that point shows the mass of the halo to which that point belongs. 

With this picture in mind, one can find analytic formulae for halo abundance as a function of their mass: the fraction 
of volume that belongs to halos in the mass range (S, S + dS) is the fraction of random walks that, starting from 
(S = 0, S = 0), first cross the barrier of height S c (z) between (S, S + dS). This can be worked out analytically (Bond 
et al. 1991) and the result is: 

f(S\S c (z))dS = -L= 6 -Mex P (- S 4f) dS (A3) 
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The total mass of the halos in this mass range and in a unit volume will then be pf(S\S c (z)). Finally, the number 
density of halos is this total mass divided by the mass of an individual halo for the S corresponding to this mass in 
the S — M relationship (See Figure A6). This gives the famous Press-Schechter formula. 

Thinking in terms of random walks has the invaluable advantage of being extendable beyond simple Press-Schechter 
formulae. Sheth & Tormen (2002) used a more realistic model of halos in which they are ellipsoidal objects instead of 
spherical, and argued that this leads to a barrier that is a function of S for a given redshift. The problem then reduces 
to a first crossing problem for a moving barrier. The result is in much better agreement with simulations. 

Lacey & Cole (1993) argue that changing the origin of the random walk from (0,0) to (Sf,S c (zf)) corresponds to 
considering only the particles that are inside a halo with mass Sf at redshift Zf and follow their history back in time. 
Then one can find the mass fraction of halos of mass Sf at redshift Zt that were part of halos within the mass range 
(Si, Si + dSi) at an earlier redshift z%. This can be achieved by simply changing the origin of the coordinate in equation 
A3 from (0,0) to (Sf, S c (zf)) (following Lacey & Cole (1993) we show 8 c (zk) as ujk)'- 



ftot(Si\uJi, Sf,w f ) dS, 
' (27r)V2(5 4 -5/) 3 / 2 



exp 



(ujj - UJf) 2 

'm-Sf) 



dSi 



(A4) 



Putting uji — uif + du> and expanding to first order in du> gives the probability that the first crossing has a jump 
from Sf to Si while one changes the height of the barrier by a tiny dui from uif to uif + dui: 

f(S f -> S l ;ui)dS l dui = ^1/2(3. _ Sf y/2 dSidiJ ( A5 ) 

Since the probability of multiple mergers is negligible for small dui, the parent halo can split only into two progenitors 
during this time interval. 
For more details see Zentner (2007) 

NUMERICAL SOLUTION 

The goal of this section is to discretize and solve Eq. (5) for f acc {Si\ui, Sf, uif) numerically using matrix methods. 
To start, we rewrite the integral equation for clarity: 

ftot(Si\uJi, Sf,Ulf) = facc{Si\0Ji, Sf,LUf) + 
pOJi f-uJi pSi — A' (Si) rSi 

/ dull I duj 2 dSl dS 2 facc(Sl\Ldl,Sf,UJf)f(S 1 —>S2;0Jl)8(uJ 1 -Ld2)ftot{Si\Ldi 7 S2,UJ2) 

Jui s Juif JS f JS 1 +A(S 1 ) 

(Bl) 

where f acc is unknown, and f merge and ftot are given in equations A5 and A4 respectively. Also, for a reason which 
will be clear later we insert an integration of the Dirac delta function into the original equation 5. Since we will 
compare this analytic result with a Monte-Carlo simulation we take (Sf,uif) as given constants which correspond to 
the mass and redshift of the parent halo. 

To solve this equation numerically wc discretize S € [Sf,Si] into N$ and ui e [uif,uii] into segments. To make 
sure that discretization is small enough to resolve mergers we demand: 

« AS (B2) 

where, like before, AS is defined as the minimum jump in S to have a merge, which for the case we consider here is 
AS = S(M(S) — M res ) — S. Also we must make sure that Aw is small enough for the probability of having more than 
one merger in that time interval be negligible. For that to be correct we demand: 

Ul ~ UJf < VAS (B3) 



We will show later that after satisfying these conditions the solution converges by making AS and Aw smaller and 
observing that the solution for the integral equation stays the same. 
We give a collective index j to any pair (Sk — kAS^ivi — IAuj): 

j = (k - 1) x N s + I (B4) 

to go back to original indices k and I we use: 



k = [ mod (j -!)] + ! (B5) 



and I can be found from this and equation B4 
Then we define these matrices: 



jijmerge _ j f merge {^ki : ; Skj , ui 3 ) if M(Sk, ) > M(Sk 3 ) + M res and uji. = uji 3 
*>} ~ \ otherwise 



(B6) 



n/r tot _ J ftot(Ski\wii, Sk^wij) if S^ > Sk 3 and LOi i > tu^ , . 

1V1 i>j ~ \ otherwise {al} 

where i,j e [1, N s x N u ]. 
Also we define this vector: 

Vl ot = f tot {S ki \uJ h ,S f ,u f ) (B8) 
All the above matrices and vectors can be computed numerically. Finally we define the unknown vector: 

Vr c = facc(SkM t ,S f ,uj f ) (B9) 
Using these definitions we can write the integral equation Bl in its discretized form as a matrix equation for vector 

ytot = yacc + M tot _ j^rnerge _ yaee ^ ,^2 ^ 



Notice that the limit of integration is taken into account in the definition of the matrices. Also notice that there is 



a factor of Aw instead of (Aw) 2 in the above formula. That's because there was a delta function in the definition of 
f merge which gives a factor of 1/Aw in discretization. Equation BIO has the solution: 

yacc = (l + M tot ■ M merge (AS) 2 Aw) ~* • V tot (Bll) 

Having found V acc , we can easily calculate the propagators involving exactly one merger, two mergers, . . . : 

yVmerge _ (AS) 2 AujV aCC • M merge ■ \F acc 
y2merge _ (/S^C/j 2 /\ ti ,\/ lmer 9 e . ^/["lerge _ yacc 

\ (B12) 

This completes our numerical solution to the integral equation for the general functions f to t and f acc . For the special 
case of spherical collapse (equations A4 and A5) a further simplification is possible by noticing that the first integration 
over S2 can be done analytically. 

MONTE-CARLO SIMULATION 

First we briefly describe how the binary merger with accretion works. For a small change in Aw the probability of 
absorbing a mass AS in this time interval is given by: 



P(AS, Aw) dS = - J )1/2 (A ^ 3/2 exp 



(Aw) 5 



2AS 



(CI) 



Starting from a parent halo with mass S p at time w we go backward in time tow - Aw. Then, if M(AS) < M res we 
consider that to be accreted mass, stop tracking its history, and take M p — M(AS) as the new parent halo at redshift 
w — Aw assuming that this new M p is larger than M res . If not, we consider that as accreted mass and do not continue 
to track its history. We choose Aw small enough to ensure that the probability of having a merger in this time interval 
becomes small. In other words we demand: 



Aw « ^S(M P - M res ) - S(M P ) (C2) 

Then we generate a random number AS consistent with the distribution CI. This is a very easy task to do since 
equation CI can be converted to a Gaussian distribution by a change of variable x = Aw/(2\/ AS*). This procedure 
is then repeated with the new halos as the parent halos until we reach the time w, where we want to compare our 
numerical results with the Monte-Carlo simulation. While making the merger tree we keep track of each halo to know 
how many times in their history they experienced a merger so we will be able to find the distribution of halos with no 
merger, one merger and so forth. 



